Determining topological order from a local ground-state correlation function 
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Topological insulators are physically distinguishable from normal insulators only near edges and 
defects, while in the bulk there is no clear signature to their topological order. In this work we 
show that the Z2 index of topological insulators and the Z index of the integer quantum Hall effect 
manifest themselves locally. We do so by providing an algorithm for determining these indices from 
a local equal time ground-state correlation function at any convenient boundary conditions. Our 
procedure is unaffected by the presence of disorder and can be naturally generalized to include weak 
interactions. The locality of these topological indices implies bulk-edge correspondence theorem. 
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I. INTRODUCTION 

The theoretical proposal and experimental discovery 
of topological band insulators pQ (TI) has been raising 
increasing interest in the condensed matter physics com- 
munity. These materials form a novel topological state 
of matter, which does not fall into the standard clas- 
sification of broken symmetries. Instead, what distin- 
guishes this phase from a trivial band insulator (BI) is 
the existence of a nontrivial Z2 topological index asso- 
ciated with the full band structure. This distinction is 
somewhat analogous to the integer quantum Hall effect 
(IQHE), whose distinguished states can be attributed to 
a Z topological index associated with full Landau levels 
of noninteracting electrons. 

Broken symmetry phases are characterized by lo- 
cal order parameters, in what is known as the Lan- 
dau paradigm [2 . It is commonly accepted that this 
paradigm does not apply to the above topological phases. 
Instead, such phases are described by the more elusive 
quantity, known as a topological order [3]. 

In IQHE the quantized Hall conductance is a direct 
manifestation of the topological order, and for the case 
of free electrons, this response function is a local bulk 
quantity [4 . As far as we know, for TI the implication of 
the topological order is only through edge effects (see, for 
example, Refs. [5HZ]) or defect-related effects [8 . There is 
no local response function that is known to characterizes 
this Z2 phase, let alone an order parameter. 

In this work we show that both the Z2 index of the TI 
and the Z index of the IQHE, as well as any gapped in- 
sulator with non-zero Chern number can be extracted 
from a local equal time ground-state correlation function. 
This implies that in these systems the topological order 
is in fact a local ground state property. It also proves 
a bulk-edge correspondence theorem for such local topo- 
logical orders. 

When interactions are taken into account, the Z index 
remains well defined [9], while it is yet unclear whether 
the Z2 index does [lOj [11]. Our formulation, however, 
remains well defined at least for weak interactions, and 
thus suggests an extension to the definition of the TI to 
weakly interacting systems. Beyond weak interactions, 



we show that either the energy gap or some "occupation 
gap" must be closed during the transition from a TI to 
a BI. 

The theoretical procedure that we provide can be 
straightforwardly adapted to form an algorithm for nu- 
merically determining the Z2 and Z indices. This algo- 
rithm is rather efficient, since one only needs to diago- 
nalize several matrices with dimensions of the order of 
the correlation length squared (2D) or cubed (3D). Such 
an algorithm may help in numerically testing a predicted 
topological order. 

II. FINDING LOCALLY THE TOPOLOGICAL 
ORDER 

Our procedure of extracting the topological index deals 
with noninteracting lattice-based band insulators. The 
required input is the equal time ground-state two-point 
correlation function 

Pij = (gs\4^\gs), (1) 

where \gs) is the many body ground state, and ipj (^) 
is the creation (annihilation) operator of an electron in 
site i (for brevity we let i encompass also the orbital 
and spin degrees of freedom). For noninteracting elec- 
trons the two-point correlation function coincides with 
the single-particle spectral projector 

Pa = E <^>Wj>> ( 2 ) 

E n <fi 

where \i) is the single-particle wave function associated 
with t/jJ, \n) is an eigenstate of the single-particle Hamil- 
tonian with an energy E n , and \i is the chemical poten- 
tial. Since we are interested in local bulk properties, we 
assume that the sites i and j reside within the bulk. 

In the following (see Property I below) we show that 
Pij is exponentially local, namely it decays exponentially 
with the distance between i and j within some correla- 
tion length l p and has only an exponentially small depen- 
dence on details of the Hamiltonian outside a local region 
around i and j. This locality allows us to discuss P^ 
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within some given A, without concerning ourselves with 
details of the edges of the region, boundary conditions, 
or the Hamiltonian outside A. In particular, for any ge- 
ometry of ^4, we can always consider a subregion with 
a geometry of a Corbino disk, or a thick torus (Corbino 
donut) in 3D, which we assume to have circumferences 
larger than l p . 

Given P on such a region, we multiply Pij by a geo- 
metric phase factor, in a way that we call "artificial flux 
insertion" 

P ij (<P) = P ij e id ^, (3) 

where </> £ R and 6ij G (— 7r,7r] is the azimuthal angle 
from site i to site j. In the proceeding we assume for 
convenience that in 2D A is of the topologically equiva- 
lent cylindrical geometry. 

A key analytical result in this work (Property II) is 
that Eq. ([3| captures the effect of real magnetic flux in- 
sertion, up to 0(e~ L / lp ) corrections, where L is the inner 
circumference of the Corbino disk (or donut). Hence the 
inclusion of flux is merely a transformation of P, which 
yields no extra information in addition to the information 
already contained in P [T2] . 

The next step is to construct ID Wannier functions out 
of P. Let X be the position operator along the open co- 
ordinate x. We define the ID Wannier functions \w n ((/))) 
to be the eigenstates of the projected position operator 
in a given flux [13] 

PXP\^\w n ((f>)) = X n ((f>)\w n ((f>)). (4) 

This definition of the Wannier functions has several ad- 
vantages. First, it relies on P rather than the Bloch wave 
functions, hence the eigenvalues are unaffected by the 
phase freedom in Fourier space, on the one hand, and it 
remains defined in the presence of disorder, on the other 
hand. Second, since P((j)) is a continuous function of 0, 
the s are also continuous in <j>. Third, in time-reversal- 
preserving systems PXP\§ is a time-reversal-invariant 
operator for <ft — 0, 7r, and thus Kramers' theorem as- 
sures that the Wannier functions come in time-reversal 
pairs with doubly degenerate eigenvalues at these points. 
Last, we prove below (Property III) that within the bulk 
each Wannier function is exponentially localized around 
its eigenvalue along x , even for cases of nontrivial Chern 
number. This implies that the x n 's of the Wannier func- 
tions that reside within the middle of the region A are 
unaffected by details of the system out of A or by the 
edges of A. 

So far we have shown that given P((j) = 0) in a local 
region A, one has sufficient data to extract x n (<j)). In 
order to extract the Z 2 index out of the x n (0)'s, we follow 
Ref. [14] and consider the pairs of eigenvalues at (j) = 
and (j) = 7r. As mentioned before, at these fluxes the x n 's 
come in degenerate pairs. The difference between a BI 
and a TI is whether the pairs at <p = are the same as 
those at (j) = tt (BI) or not (TI), as depicted in Fig. [I] 

The Z index is even simpler to extract, since there 
is no degeneracy in the x n 's, and all the ID Wannier 



functions move in the same direction [15]. According to 
gauge invar iance, the eigenvalues at <\> = are the same 
as those at (j) = 2ir, but each x n may be carried with the 
flux to x m . If the labeling of the eigenvalues is such that 
x n+ i > x n , then the Z index equals to ra — n, as depicted 
in Fig. [2] 

The four Z 2 indices of the 3D TI {y^v x ^y^ v z) can 
be extracted by generalizing the pair switching criterion 
[IS] . Given a sample with periodic boundary conditions 
in y and z, the Wannier centers are carried with two inde- 
pendent fluxes x n (cj)y, cj) z ), where <\> y (cj) z ) corresponds to 
the phase twist in y (z). Now v y = 1 if the pairs switch 
partners between ((j) y ,(j)z) = C^O) and (tt, tt) and, ac- 
cordingly, v z = 1 for switching between (0,7r) and (71", tt). 
If the pair switching in the course (0, 0) — » (0, tt) differs 
from that in (tt, 0) —> (tt : 7r), then = 1. In order to 
extract v x we must have P also in a geometry that is 
periodic in x and y. 

We can see that in 3D TI the geometry on which P is 
given plays an important role. The 3D generalization of 
the 2D cylinder is periodic boundary conditions in two di- 
rections, while the 3D generalization of the Corbino disk 
is the Corbino donut. The Corbino donut can be iso- 
lated from any given sample by discarding sites from the 
projector, and therefore properties inferred from these 
geometries must be local and isotropic. On the contrary, 
the 3D cylinder cannot be isolated form larger samples, 
and therefore properties inferred from this geometry need 
not be local nor isotropic. In the Corbino donut we con- 
sider only the flux that resides inside the donut, and X 
denotes the open radial coordinate. For vq = 1 we expect 
pair switching, since gapless states will appear at the sur- 
face [16]. On the other hand, since a Corbino donut can 
be isolated from a larger sample in any orientation, it is 
impossible to extract from it the orientation-dependent 
weak indices. 

In the presence of interactions, the ground- state two- 
point correlation function remains well defined. However, 
it no longer corresponds to a projection matrix, which is 
the requested input to the process. Prior to the addi- 
tion of interactions all the eigenvalues of the correlation 
function were either or 1, corresponding to occupied 
and unoccupied states. One can think of this as an oc- 
cupation gap of exactly 1. As interactions are gradually 
increased, and as long as there is no phase transition, the 
correlation function is expected to change continuously, 
and the occupation gap to remain finite. Provided that, 
it is possible to extract a truncated projector out of the 
correlation function in a way which is independent of the 
truncation process, as will be proven below. 

In general, locality of topological order implies a bulk- 
edge correspondence theorem, provided that the topolog- 
ical order is well defined for insulating phases. Assume 
that two samples which belong to two different classes 
are attached. If the entire system remains gapped, it be- 
longs to a single topological class. However, according 
to the locality of the order, the system seems to belong 
to two classes simultaneously. This contradiction implies 
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FIG. 1: The centers of the ID Wannier functions Xn cLS a func- 
tion of the artificial flux (j) of the 2D Kane and Mele model, 
as extracted from the local two-point correlation function at 
(j) = 0, of both a topological (left) and a trivial (right) insu- 
lators, (a) In the TI the centers switch pairs between = 
and cj) — 7r. (b) In the BI the centers are essentially fixed, 
(c) Zooming-in shows that they actually move, but without 
switching pairs. 



a closure of the gap, which can only take place at the 
interface. In the other way, if the gap can remain open 
between the two attached samples, the classifying order 
must be nonlocal. For example, weak topological insula- 
tors may have a gapped surface (in a stacking direction) 
[16], which means that the weak order is nonlocal, in 
agreement to what has been stated above. 



III. NUMERICAL IMPLEMENTATIONS 

The above approach can be easily adapted to a com- 
puter algorithm. For a given two-point correlation func- 
tion P, one can diagonalize the matrix PXP\§ for var- 
ious fluxes, and examine the resulting eigenvalue spec- 
trum. However, due to finite size effect, some of the 
eigenstates of PXP would be localized near the edges 
of the region A. These eigenstates are edge dependent 
and, therefore, do not reflect the physical behavior of the 
bulk. Fortunately, due to the exponential localization 
of the eigenstates, one can easily distinguish these state 
from the bulk eigenstates and discard them, as proven in 
the Appendix. 

We have carried out this process on the 2D Kane and 
Mele model of the TI [17] in the presence of disorder. 
The Hamiltonian of this model is parameterized by near- 
est neighbor hopping t, staggered on-site potential \ V1 S z 




FIG. 2: The centers x n versus the artificial flux of a single 
spin copy of the quantum spin Hall effect model, as extracted 
from P at zero flux. Each center replaces its proceeding center 
while carried with the flux, x n (27r) = x n +i(0), in all scales, 
as expected from v — 1. 



conserving spin-orbit interaction \so, and Rashba spin- 
orbit interaction A#. We took a cylindrical 18x42 lattice, 
with t 1,A 50 0.1, 0.05, for the BI X v = 0.9, 

while for the TI X v = 0.1. Both Hamiltonians were sub- 
ject to a uniformly distributed random potential of mag- 
nitude 0.1. The spectral flows x n {(j)) at the middle of the 
cylinder of both the trivial and topological phases are 
depicted in Fig. [T] The difference in the pair switching 
behavior is clearly visible. 

Similarly, we performed the procedure on a single spin 
copy of the quantum spin Hall effect model [18 , which is 
equivalent to IQHE. The parameters of this models are 
the hopping element of the two bands B — D and B + D, 
the interband hopping A, and the energy gap M. Figure 
[2] depicts the movement of the ID Wannier functions at 
the middle of a cylindrical 30 x 40 lattice, with 5 = 1, 
D = 0.2, A = 0.4, and M 0.1, which yield v = 1. A 
uniform disorder of magnitude 0.3M is also present. It 
is evident that x n (27r) = x n+ i(0), which means that this 
insulator belongs to class 1 of the Z index, as expected. 

A comment is in order that different x n (0)'s may ap- 
pear to cross at accidental values between <\> — and 
(j) = 7r; for example, see Fig. [ljc). In such cases one 
should relate the upper state before the crossing to the 
upper state after the crossing and the same for the lower. 
This process is equivalent to opening a gap at the crossing 
point by some small local perturbation. Moreover, due 
to Winger noncrossing theorem [19] , this gap is probably 
there, only it is not visible within the numerical accuracy. 



IV. PROOFS OF ANALYTICAL PROPERTIES 

Now we turn to prove the three properties that have 
been stated above. The band insulator lattice Hamilto- 
nian H is characterized by a spectrum with a finite gap 
A around some value /i, and a maximal energy in ab- 
solute value D above fi. We also assume that H does 
not couple two sites which are more than 1^ sites apart. 
Since we are interested in bulk properties, we assume that 
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the periodic dimensions of the lattice are of finite size L, 
while the open coordinate is infinite (finite size systems 
are discussed in the Appendix). 

The starting point of the proofs is to develop a repre- 
sentation of the projector P as a finite polynomial in H. 
The projector P can also be expressed as 



lim - 



1 - erf 



H 



eA 



(5) 



where erf(x) denotes the Taylor series of the error func- 
tion. The equivalence of this expression to the definition 
can be easily verified in the eigenbasis of H. Accordingly, 
we define the approximate projector P e to be the same as 
P but with finite e«l. If the error of the approximation 
is measured by the Euclidean norm, we can bound it by 
||P - P e || = (1 - erf(e-,)/2 < e^ 2 . 

Since the erf(x) is an entire function, P e can be ap- 
proximated by taking only the N first terms of the series 



N,e 



-\i-—y 



^ ^ n!(2n + l) V eA 

v n=0 v ' x 



(6) 



According to Taylor's theorem there exists xq G 
(0, (D/A)e~ 1 ) such that \\P 6 -PnA = ^ 2iV+3 /[vW + 
1)!(27V + 3)]. By using Stirling's approximation, and 
choosing N e = e 2 (D / A) 2 e -2 ^> e -2 , we can bound 



Pe-P N J < 



D ( eD 



2 \ N+l 



NAe \NA 2 e 2 



(7) 



Consequently, we can conclude that P/v e , e is an excel- 
lent approximation of P, with an error of ||P — P/v e , e || < 
e -1 / e = e _< ^, with Q as the measure of the accuracy. 
Therefore, in order to approximate P with accuracy Q, 
it is sufficient to take Nq = e 2 (D/ A) 2 Q terms. 

Property I: P^ decays exponentially with the dis- 
tance between i and j within some correlation length l pi 
and has only an exponentially small dependence on the 
Hamiltonian outside a local region around i and j. 

Proof: According to the assumptions, [H 2N+1 )ij van- 
ishes for two sites i,j which are separated a distance 
Tij > (2N + l)lh ~ 2Nlh. Accordingly, [P/v, e ]ij also van- 
ishes, while Pij may be of order e - ^^. Therefore, for 
a given > 2/^, we choose N = [r^/2^], which keeps 
[PN,e]ij as zero, while giving a maximal Q(N). Following 
this choice one finds that 



\Pij\ <e- r »/ 1 *, 
l p = 2e 2 (D/A) 2 l h . 



(8) 
(9) 



Moreover, P^ for two sites within a ceratin region A 
does not depend on the matrix elements H^i of two sites 
within region B as long as the distance between A and B 
is much larger than l p . This can be shown by choosing 
N = rjuslllh, where rjsj$ is the minimal distance between 



A and B. Now [P/v, e ]ij does not depend on Hm, and P^ 
may depend on it on the order of e~ TAB ^ lp [20] . 

This way of creating an exponentially local projector 
from the Hamiltonian can be performed on any gapped 
matrix with a finite band width. In particular the two- 
point correlation function of a weakly interacting system 
can be deformed into a projector, as required for the 
index determining procedure, as long as the occupation 
gap remains finite. 

Property II: Magnetic flux which threads a cylinder 
or a torus affects P^ by a simple geometric phase factor. 

Proof: Insertion of the magnetic flux <I> affects H in 
the uniform gauge by 



H H (<S>) = H kl e i2 ^' L *\ 



(10) 



where y k \ G (— L/2,L/2] is the azimuthal distance from 
site k to site and and $0 is the flux quanta. Consider 
Hki within a region C that cover less than a half of the 
circumference of the cylinder or the torus. One can adopt 
a convention for the coordinate y G (0, L) that avoids the 
branch cut inside C and write y k i = y k — y\. Within C 
Eq. ( 10 ) can now be written as a unitary transformation 



H{$)\c = U{$)HU\$)\c> 
U kl ($) = 5 kl e i2 *y^ L *°, 



(11) 

k,leC. (12) 



Given P^-, we can approximate it by N = L/Alh terms. 
Since [P/v, e ]*j does not vanish only for r^ ~ l p « L, it 
depends on Hj k only for j, k within a region of size P/2, 
which covers not more than a half of the circumference. 
Hence we may substitute Eq. ( pTj ) in Eq. (|6| and obtain 
Eq. Q by identifying <fi = $/$o and % = 

Property III: The ID Wannier functions, defined as 
eigenstates of the operator PXP, are exponentially local- 
ized in the x direction around the corresponding eigen- 
values. 

Proof: Consider the action of the projector P on some 
normalized state \a) which is localized around x a within 
a width l a . Since P^ is local with width l pi P\a) is still 
localized at the vicinity of x a: but it might be as wide 
as I a + l p . Note that P\a) is not necessarily normalized 
to 1 but may have a smaller norm. Since the position 
operator X is diagonal in the position basis, (X — x a )\a) 
is localized almost in the same manner as |a), but its 
norm may increases up to 21 a . 

Following this, if we begin with a particle at site z, 
and apply P(X - x t )P on it, then \\P{X - Xi)P\i)\\ < 



Ml 



< 



2l p . Applying it M times yields ||(P(X - xi)P) 
(2lp) M M\ « V2^M(2l p M/e) M . 

Now let \w n ) be an eigenstate of PXP with eigen- 
value x n ^ 0, and consider \i) with \x n — Xi\ > 2l p . 
It can be shown that P\w n ) = \w n ), yielding (i\w n ) = 
(x n - Xi)- M (i\(P(X - Xi)P) M \wn). Following the in- 
equality \(a\b)\ < |||a)|| • |||6)||, we have \(i\w n )\ < 
V / 27rM [2l p M/(e(x n — Xi))] M . For given x n and Xi we 
choose M = [\xi — x n \/2l p ], and get the exponential lo- 
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calization of the ID Wannier function [21] 

|<*K>| < <J*\xi-x n \ll p e-\ x *- x »U 21 *. (13) 



V. CONCLUSION 

To conclude, in this work we proved that the Z and 
Z2 topological indices can be extracted from the ground- 
state equal-time two-point correlation function at zero 
flux, given on any section larger than some correlation 
length. This implies that the order in such topological 
phases can be thought of as a local ground-state property. 
In the case of 3D TI it was suggested that the strong 
index is indeed local in the above sense, while the weak 
indices carries some global information. 

Heuristically, one can reach a similar conclusion using 
entanglement spectrum [22]. Indeed, on a finite re- 
gion is related to the reduced density matrix. A gapless 
spectrum of Pij is therefore an indication of nontrivial 
topology [23j[23]. Nevertheless, given a finite region one 
cannot differentiate a gapped spectrum from a gapless 
one without inserting fluxes. 

It will be interesting to establish this approach in in- 
teracting systems and to apply it on fractional quantum 
Hall states. Preliminary numerical results indicate that 
the electron two-point correlation function of the ground 
state at filling factor 1/3 is proportional to the one of 
filling factor 1. Hence, it may be possible to extract the 
topological order of fractional states in a similar way. 

We hope that our viewpoint will encourage the efforts 
for finding manifestations of the Z2 index through local 
bulk properties. 
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Appendix A: Edge effects 



The starting point of the proofs in Sec. [TV] was expand- 
ing the spectral projector P as a series in powers of the 
Hamiltonian iJ, with the spectral gap A as a control pa- 
rameter. An edge may give rise to gapless edge states 
and so our first task is to establish the expansion in the 
presence of such states. 

We assume that H is characterized as before, only now 
the spectrum of H includes both bulk states with energy 
greater than the gap, and edge states with subgap energy, 
which are exponentially localized along the edge. The 
projector P can than be divided into bulk and edge parts 



P = P 



bulk 



(Al) 



]T \n)(n\+ l n >H> 

E n <fi-A fi-A<E n <fi 



where \n) denotes an eigenstate with eigenvalue E n , and 
fi is the chemical potential, as before. Since P ed 9 e is 
composed only of the edge states, it decays exponentially 
into the bulk. Therefore [P edge ]ij is exponentially small, 
if either i or j (or both) reside within the bulk. 

In a similar manner to what we have done above, we 
introduce an error function approximation to the projec- 
tor 



p pbulk _|_ pedge 



(A2) 



?J('- rf ( 



£„<M-A 



eA 



\n)(n\ 



fi-A<E n <fi 



E n ~ H 

eA 



\n)(n\. (A3) 



We can see that IIP' 



>bulk 



ybulk I 



< e 



-1/e 2 



as before, 



while P^ d 9 e is a poor approximation to P eu ^° 5 since 
the summation coefficients spread form to 1. Indeed 
pedge j g a p 00r approximation of P ed 9 e at the edge. 
However, within the bulk both [Pf ge ]ij and [P ed9 %j 
are exponentially small, and therefore also the error 

^pedge _ pedge^.y 

This implies that for bulk-bulk or bulk-edge correla- 
tion, P may still be approximated by P e within an expo- 
nential accuracy. P e can be Taylor expanded up to some 
finite order N, and by choosing optimal values for e and 
N we can bound the error. Since this part of the proof 
is unaffected by the edges, we refer the reader back to 



Sec. [TV] We turn to discuss the effect of edges on prop- 
erties I— III. 



In the main text we avoided edge effects by consider- 
ing geometries with periodic dimensions and an infinite 
open coordinate. In practice, a finite open coordinate 
must be used, causing numerical artifacts that should be 
filtered out in order to reveal the bulk behavior. This 
can be done in a controlled manner provided that the 
ID Wannier functions of the bulk and the edge are dis- 
tinguishable. The first part of this Appendix establishes 
this distinction, and the second part presents the actual 
numerical procedure. 



Property I (the exponential locality of P) relies on the 
serial expansion, and therefore is valid only within the 
bulk. Nevertheless, since the edge states decay into the 
bulk, P decays exponentially also at the vicinity of the 
edge but only perpendicularly to the edge. 

Property II (the artificial flux insertion) is valid as long 
as P decays along the periodic coordinates, which does 
not necessarily happen in the presence of edge states. 
Thus the artificial flux insertion approximation is valid 
only far from the edge. 
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Property III states that the ID Wannier functions are 
exponentially localized around their eigenvalues along the 
open coordinate X. As long as the ID Wannier functions 
are produced by diagonalizing PXP, and P is a true pro- 
jection operator over the entire physical system, the proof 
remains unchanged in the presence of edges. Neverthe- 
less, we do not wish to limit ourselves to cases in which 
the entire P matrix is known, since this is not a local 
quantity. Our algorithm uses P which is given on some 
local area with a geometry of a cylinder or a Corbino 
disk. But P which is truncated to some local region is 
generally not a projection matrix. 

The projection property P 2 = P was used in the orig- 
inal proof only once, when it was stated that it can be 
shown that P\w n ) = \w n ), where \w n ) is an eigenfunc- 
tion of PXP with eigenvalue x n ^ 0. This property was 
required in order to establish the equality 



<*K> = (x n - x t )- M (i\(P(X - Xi )P) M K>, (A4) 

where \i) is a state localized at site i . If the truncated 
P is no longer a projector, then [PXP, P] ^ 0, and \w n ) 
may not be an eigenstate of P. But since we are inter- 
ested in bulk properties, it suffices to prove that Eq. ( A4) 



is valid far from the edges of P (which may differ from 
the physical edges due to the truncation). 

It is useful to divide the the system into three groups: 
S will denote the entire system, A C S the truncation 
area, and B C A the inner region of A, which will be 
defined by the set of all sites which are far form the edges 
of A, on the scale of the the correlation length l p . Recall 
that P is a piece of the true projector of the entire system, 
which we denote by P. Following the locality of P, we 
find that for i or j in B 



\P ]ij — ^ ^ Pik-Pkj — ^ ^ Pik-Pkj 
keA keA 



^ PikP} 



kj — *ij 



(A5) 



kes 



up to corrections that are exponentially small in the dis- 
tance between the site in B and the edge of A : divided by 
l p . We can see that P is still a projector up to boundary 
effects. 

Since P is approximately a projector within 6, 



(i\P\w n ) = (l/x n )(i\PPXP\w n ) 
*(l/x n ){i\PXP\w n ) = 



(A6) 
Vi e B, 



for x n ^ 0. This means that as far as region B is con- 
sidered, \w n ) is indeed an eigenstate of P. Seemingly 
we recovered Eq. ( |A4[ ) for any site i within B. However, 
(P(X — Xi)P) M \i) may b e as wide as M • 2Z P , which re- 
stricts the validity of Eq. ( |A4| ) to M < d(i, B)/2l p , where 
d(i,B) is the distance between site i and the edge of B. 

The exponential localization of \w n ) around x n was 
achieved by choosing M = [\xi— x n \/2l p ]. The restriction 
on M is then translated to a restriction on the exponen- 
tial localization to be valid only for \x{ — x n \ < d(i,B), 
although Xi and Xji are both in B. 



So far we have shown that for x n in B, its eigenfunc- 
tion \w n ) decays exponentially with the distance. Note 
that this does not exclude the possibilities that \w n ) re- 
sides at the edge of ^4, or both at the edge and around x n . 
The last scenario is, however, highly nongeneric. Assume 
that some \w n ) is indeed doubly localized in that fashion. 
Up to exponential accuracy we may split it into sum of 
two functions \w n ) = \w n ,B) + \w n ,A), where \w ni B) is 
localized around x n and \w n ,A) is localized around the 
edge. Since PXP is local, it does not couple these two 
functions, and \w n ) = \w n ,B) — \w ni A) is also an eigen- 
function, with an eigenvalue that is exponentially close to 
x n . Such an almost degeneracy is of course nongeneric. 





FIG. 3: Demonstration of bulk-edge separation in the single 
spin copy of the quantum spin Hall effect. The projector P 
was created by truncating a 30 x 34 cylinder, denoted by A, 
from the 30 x 40 torus. The full spectrum of PXP\ ct) (top, 
black) is trivial both in the bulk and at the edge. Omitting 
states that are localized at the edges gives a pure bulk spec- 
trum at region B (bottom, red), which is indeed nontrivial. 



To conclude, we have seen that edges do not affect 
the bulk properties of P. Moreover, the ID Wannier 
function with eigenvalue in the bulk are localized either 
around their eigenvalue or at the edge. Therefore, given 
the correlation function P on an arbitrary geometry, one 
can isolate/truncate a cylinder or a Corbino disk from 
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it and distinguish the ID Wannier functions that reside 
within the bulk . Due to their localization properties, the 
ID Wannier functions are unaffected by the edges or the 
truncation process, and the motion of their eigenvalues as 
a function of the flux is therefore purely a bulk property. 

We now demonstrate how these analytical statements 
are exploited in the computer algorithm. For that pur- 
pose we fully diagonalized the Hamiltonian of the single 
spin copy of the quantum spin Hall effect, for a periodic 
30 x 40 lattice, with the parameters that are given in 
Sec. |III| The projector P was than produced by sum- 
ming all projectors on the states with negative energy 
and was characterized by a correlation length of approx- 
imately 2. Region A was chosen be a 30 x 34 cylinder 
out of the full torus, and all the elements P^ with i and 
j out of A were omitted. 

The next step was to construct the spectrum {x n ((/>)} 



by diagonalizing PXP|^, where X in region A ranged 
from 4 to 37. The full spectrum is trivial, as depicted at 
the top of Fig. [3] At the edges of A edge states hybridize 
with the outmost bulk states, and gaps are opened. At 
the center of A most eigenvalues belong to bulk states, 
while the branch of eigenvalues that crosses the spectrum 
from side to side belongs to an edge state. 

In order to remain with bulk effects only, we used 
the localization property of the ID Wannier functions, 
which guarantees that the bulk wave function are local- 
ized within the middle of region A, denoted above by B. 
We chose B to be the central cylinder of 20 sites, and 
discarded all the ID Wannier function that more than 
5% of their weight is outside B. In this way we assured 
that all the states in B are purely bulk states, and the 
nontrivial nature of the spectrum became apparent, as 
seen in the bottom of Fig. [3] 
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